Introduction
The purpose of this document is to provide a reproducible record of all analyses and figures in the main article. The main article is focused on the effect of response diversity on community stability in fluctuating environments. We are going to look at the effect of response diversity, richness, temperature and nutrients on community temporal stability. Specifically, we are going to look at the effect of fundamental balance (our measurement of stability) on temporal stability. Then, as response diversity is thought to stabilize temporal stability of aggregate community properties via asynchrony, we are going to look at the relationship between response diversity and asynchrony.
Finally, as multiple evidence suggests that compensatory dynamics and temporal stability are determine by species interactions, we are going to analyse the effect of species interactions on stability to understand if they are more important than response diversity in driving temporal stability of community biomass.
This document is produced by an Rmarkdown file that includes code to reproduce from data all results presented in the main article.
Load datasets, Data wrangling and balance calculation
divergence_df <- read_csv("Data/divergence_df.csv")
load("Data/dens_biomass_poly.RData")
dd_all_pred<-read.csv("Data/morph_dd_pred.csv")
dd_all_pred_nonoise<-read.csv("Data/morph_dd_pred_nonoise.csv")
load("Data/ciliate_traits.Rdata")
df_slopes <- read_csv("Data/df_slopes_cor.csv")
# needs to have id_new variable
ciliate_traits <- ciliate_traits %>%
dplyr::mutate(
# Remove dots from the date
cleaned_date = gsub("\\.", "", date),
# Extract the part of id after the underscore
id_suffix = sub(".*_(.*)", "\\1", id),
# Combine cleaned_date, id_suffix, and species_initial into a new variable
id_new = paste0(cleaned_date, id_suffix, composition)
) %>%
# Optionally, remove the intermediate columns to clean up
dplyr::select(-cleaned_date, -id_suffix,-new_id)
uniqueN(ciliate_traits$id_new)==nrow(ciliate_traits) # all unique ;)
id_dd<-full_join(dd_all_pred,dplyr::select(ciliate_traits,id_new,biomass),join_by("id_new"))
## add day variable
#create a day variable from the date variable
id_dd<-dplyr::mutate(id_dd,date=as.Date(date,format = "%d.%m.%y"))
earliest_date<-min(id_dd$date)
days_since_earliest<-as.numeric(id_dd$date-earliest_date)+1
id_dd<-id_dd%>%dplyr::mutate(day=days_since_earliest)
#create a summarised df on microcosm level with each species seperate
# Make sure, that we have n_frames and not N_frames
names(id_dd)[names(id_dd) == "N_frames"] <- "n_frames"
#extrapolation_factor <- 9.301902 # for 16 x magnification
extrapolation_factor <- 9.828125 # for 25 x magnification
video_biomass_species <- c( "C", "P", "S","D","L","T")
biomasses <- id_dd %>%
dplyr::group_by( day,temperature,nutrients,sample_ID,composition,predict_spec) %>% # group by xxx
dplyr::summarize(
biomass = sum(biomass * n_frames, na.rm = TRUE) / (1 * 125) # if not 3 videos corrections is done below with dens_factor
) %>%
dplyr::mutate(
biomass = biomass * extrapolation_factor,
)
biomasses<-biomasses%>%dplyr::mutate(biomass=biomass*1000)
dd_ts_id<-biomasses
#fill up missing dates with biomass<-0
fill_dd<-expand.grid(sample_ID=unique(dd_ts_id$sample_ID),day=unique(dd_ts_id$day),predict_spec=unique(dd_ts_id$predict_spec))
complete_ts<-full_join(fill_dd,dd_ts_id,join_by(sample_ID,day,predict_spec))
complete_ts$biomass[is.na(complete_ts$biomass)]<-0
complete_ts<-complete_ts%>%dplyr::mutate(composition=sub("_.*", "", sample_ID))
complete_ts<-complete_ts %>%
dplyr::mutate(temperature = sapply(strsplit(as.character(sample_ID), "_"), function(x) paste(x[3], x[4], sep = "-")))
complete_ts<- dplyr::mutate(complete_ts,nutrients = gsub(".*Nut(.*?)_.*", "\\1", sample_ID))
# Now remove wrong combinations of composition and predict_spec / predict_spec
complete_ts<- complete_ts %>%
rowwise() %>%
dplyr::filter(predict_spec %in% unlist(strsplit(composition, ""))) %>%
ungroup()
complete_ts<-dplyr::mutate(complete_ts,temperature=as.character(temperature),
nutrients=as.character(nutrients),
richness=nchar(composition))
complete_ts<-complete_ts%>%group_by(sample_ID,composition,day)%>%dplyr::mutate(tot_biomass=sum(biomass))
complete_ts<-complete_ts%>%dplyr::mutate(biom_contribution=biomass/tot_biomass)
df_biomass_mod <- complete_ts
complete_ts<-complete_ts%>%dplyr::mutate(temperature=paste0(temperature," °C"),
nutrients=paste0(nutrients," g/L"))
# introduce slopes of
names(df_slopes)[names(df_slopes)=="species_initial"]<-"predict_spec"
slope_ts<-full_join(dplyr::select(df_slopes,nutrients,predict_spec,temperature,slope),complete_ts)
slope_ts<-slope_ts%>%dplyr::mutate(w_slope=biom_contribution*slope,
sign=sign(slope))
slope_ts<-slope_ts%>%group_by(sample_ID,temperature,nutrients,richness,composition,day,tot_biomass)%>%dplyr::summarize(
sum_w_slopes=abs(sum(w_slope)),
mean_abs_slope=mean(abs(slope)),
sum_abs_slope=sum(abs(slope)),
abs_sum_slope=abs(sum(slope)),
symmetry=abs(sum(sign)))
slope_ts<-slope_ts%>%dplyr::mutate(richness=as.factor(richness))
##create new variable where it checks, where the last observation =0 is; with complete_ts
aggr_ts <- slope_ts %>%
group_by( sample_ID) %>%
arrange(day) %>%
mutate(
# Create a flag for non-zero tot_biomass
non_zero_biomass = tot_biomass != 0,
# Find the last non-zero day
last_non_zero_day = ifelse(any(non_zero_biomass), max(day[non_zero_biomass], na.rm = TRUE), NA),
# Find the first zero day after the last non-zero day
first_zero_day = ifelse(
!is.na(last_non_zero_day),
min(day[!non_zero_biomass & day > last_non_zero_day], na.rm = TRUE),
NA
),
# Flag for days after the first zero day
is_after_first_zero_day = ifelse(!is.na(first_zero_day), day > first_zero_day, FALSE)
) %>%
ungroup()
aggr_ts<-aggr_ts%>%mutate(rep_var=sub("_[^_]+$", "", sample_ID))
biomass_ts<-aggr_ts%>%group_by(day,temperature,nutrients,richness)%>%summarize(tot_biom=mean(tot_biomass),se_tot_biom=sd(tot_biomass)/sqrt(as.numeric(length(tot_biomass))))
Biomass
Let’s have a look at the biomass dynamics in the different environmental treatments.
tot biomass plot

Figure 1 : Community total biomass during the experiment in different environmental treatments. Different color represent richness levels.
Main Results
We now look at the main results of the experiment. We are going to look first at the effect of richness, temperature and nutrients on community temporal stability. Then, we are going to look at the relationship between divergence (original response diversity metric) and temporal stability. Finally, we are going to look at the relationship between response diversity and temporal stability.
In the whole analysis, we calculated the temporal stability of total community biomass as the inverse of the coefficient of variation (ICV) (i.e. \(\frac{\sigma}{\mu}\)).
Effect of T, N and R
Figure 2: Effects of richness (a), temperature (b), and nutrients (c) on community total biomass temporal stability.
We can see that richness does not have a clear effect on community temporal stability, while stability was higher at lower temperature, and nutrients increased community temporal stability.
Effect of Divergence
We look at the relationship between divergence (our original response diversity metric) and stability

Figure 3: Relationship between Divergence and temporal stability of total community biomass.
Divergence is positively related to temporal stability, suggesting that response diversity promotes stability. However, the relationship between divergence and stability becomes weaker as richness increases. We think that this is due to divergence considering only the responses of the 2 most “responding” species. Thus, when species richness increases, disregarding the responses of the other species in the community except the 2 responding the most makes the relationship between response diversity and stability weaker.
This is why, after running the experiment, we developed another metric to measure response diversity, which we called balance, and that is presented in the main text of the publication.
Balance has several desirable features that makes it a more suitable metric than divergence: Independence of richness, higher predictive power, and accounts for the responses of all species in the community (as opposed to divergence that accounts for only the 2 most “responding” species).
Here, we provide extensive evidence of why balance is a better metric to measure response diversity than divergence, and thus justifying focusing the analysis around balance.
Comparing Divergence and Balance
Predictive power of Divergence and Balance
We first compare how well divergence and balance predict stability (predictive power).
Balance
#
mod1 <- lm(data=complete_aggr,log10(stability)~log10(balance_f))
# Check model assumptions
#check_model(mod1)
Divergence
mod2 <- lm(data=complete_aggr,log10(stability)~(divergence))
# Check model assumptions
#check_model(mod2)
Table 1: Comparison of model performance of divergence and balance as predictors of stability. Model 1 has balance as predictor and model 2 has divergence as predictor.
|
model
|
AIC
|
AICc
|
BIC
|
R2
|
R2_adjusted
|
RMSE
|
Sigma
|
|
1
|
-89.27328
|
-89.17286
|
-78.79409
|
0.1917679
|
0.1884142
|
0.1510344
|
0.1516599
|
|
2
|
-55.71579
|
-55.61538
|
-45.23661
|
0.0720796
|
0.0682293
|
0.1618316
|
0.1625017
|
A model with Balance as predictor performs better than one with divergence as predictor, and it explains more of the variance in stability than divergence.
Moreover, from Figure 3, it looks like divergence declines in performance as richness increases. Let’s test this analytically.
To do than we build a linear model having stability as response variable and either log10(balance) or divergence as predictor for each richness level. We then extract the R squared of the models and their standardised estimates. (standardized estimates were calculated centering divergence and balance using the function scale()).
# getting model estimates for each richness level
lm_divergence_richness_E <- complete_aggr %>%
nest(data = -richness) %>%
mutate(
model = map(data, ~ lm(log10(stability) ~ scale(divergence), data = .x)),
results = map(model, broom::tidy)
) %>%
unnest(results) %>% dplyr::filter(term=="scale(divergence)")
# getting model R squared for each richness level
lm_divergence_richness_R <- complete_aggr %>%
nest(data = -richness) %>%
mutate(
model = map(data, ~ lm(log10(stability) ~ scale(divergence), data = .x)),
results = map(model, broom::glance)
) %>%
unnest(results)
# getting model estimatesf or each richness level
lm_balance_richness_E <- complete_aggr %>%
nest(data = -richness) %>%
mutate(
model = map(data, ~ lm(log10(stability) ~ scale(log10(balance_f)), data = .x)),
results = map(model, broom::tidy)
) %>%
unnest(results) %>% dplyr::filter(term=="scale(log10(balance_f))")
# getting model R squared for each richness level
lm_balance_richness_R <- complete_aggr %>%
nest(data = -richness) %>%
mutate(
model = map(data, ~ lm(log10(stability) ~ scale(log10(balance_f)), data = .x)),
results = map(model, broom::glance)
) %>%
unnest(results)
Figure 4: Performance comparison of divergence vs balance. In (a), the R squared of linear models for divergence and balance are shown for each richness level. In (b), the estimates of the linear models for divergence and balance are shown for each richness level.
We can see that the R squared of divergence as predictor of stability becomes smaller as richness increases, while the R squared of balance as predictor of stability does not (actually increases slightly).
Comparing unique explanatory power of balance and divergence
Now we build a linear model were stability is modeled as a function of balance and divergence.
Then, we compared the variance explained by the full model compared to a model containing either only balance or only divergence.
Full model - balance and divergence
lm_div_balance <- lm(data=complete_aggr,log10(stability)~log10(balance_f)+divergence)
# Check model assumptions
# check_model(lm_div_balance)
model with only divergence
lm_div <- lm(data=complete_aggr,log10(stability)~divergence)
# Check model assumptions
# check_model(lm_div)
model with only balance
lm_balance <- lm(data=complete_aggr,log10(stability)~log10(balance_f))
# Check model assumptions
# check_model(lm_balance)
Comparision full model vs divergence only and balance only
Table 2: Comparison of model performance of divergence, balance and both as predictors of stability. Model 1 has both balance and divergence as predictors, model 2 has divergence as predictor, and model 3 has balance as predictor.
|
model
|
AIC
|
AICc
|
BIC
|
R2
|
R2_adjusted
|
RMSE
|
Sigma
|
|
1
|
-88.97683
|
-88.80876
|
-75.00458
|
0.1974141
|
0.1907259
|
0.1505060
|
0.1514437
|
|
1
|
-55.71579
|
-55.61538
|
-45.23661
|
0.0720796
|
0.0682293
|
0.1618316
|
0.1625017
|
|
2
|
-89.27328
|
-89.17286
|
-78.79409
|
0.1917679
|
0.1884142
|
0.1510344
|
0.1516599
|
Comparision full model vs balance only
Table 3: Anova table: a model with both balance and divergence as predictors is not significantly different from a model with only balance as predictor.
anova1 <- anova(lm_div_balance, lm_balance)
# Convert to tidy format
anova_tidy1 <- broom::tidy(anova1)
# Display the tidy ANOVA table using gt with formatted p-values and adjusted size
anova_tidy1 %>%
gt() %>%
cols_label(
term = "Term",
sumsq = "Sum of Squares",
df = "DF",
statistic = "F Statistic",
p.value = "p-value"
) %>%
fmt_number(
columns = vars(p.value),
decimals = 3
) %>%
tab_style(
style = cell_text(weight = "bold"),
locations = cells_body(
columns = vars(p.value),
rows = p.value < 0.05
)
) %>%
tab_options(
table.width = px(800), # Adjust table width (e.g., 400px)
table.font.size = px(12), # Adjust font size (e.g., 12px)
data_row.padding = px(10) # Adjust row padding (e.g., 4px for more compact rows)
)
| Term |
df.residual |
rss |
DF |
Sum of Squares |
F Statistic |
p-value |
| log10(stability) ~ log10(balance_f) + divergence |
240 |
5.504447 |
NA |
NA |
NA |
NA |
| log10(stability) ~ log10(balance_f) |
241 |
5.543171 |
-1 |
-0.03872444 |
1.688429 |
0.195 |
Comparision full model vs divergence only and divergence only
Table 4: Anova table: a model with both balance and divergence as predictors is significantly better from a model with only divergence as predictor.
anova2 <- anova(lm_div_balance, lm_div)
anova_tidy2 <- broom::tidy(anova2)
# Display the tidy ANOVA table using gt with formatted p-values and adjusted size
anova_tidy2 %>%
gt() %>%
cols_label(
term = "Term",
sumsq = "Sum of Squares",
df = "DF",
statistic = "F Statistic",
p.value = "p-value"
) %>%
fmt_number(
columns = vars(p.value),
decimals = 3
) %>%
tab_style(
style = cell_text(weight = "bold"),
locations = cells_body(
columns = vars(p.value),
rows = p.value < 0.05
)
) %>%
tab_options(
table.width = px(800), # Adjust table width (e.g., 400px)
table.font.size = px(12), # Adjust font size (e.g., 12px)
data_row.padding = px(10) # Adjust row padding (e.g., 4px for more compact rows)
)
| Term |
df.residual |
rss |
DF |
Sum of Squares |
F Statistic |
p-value |
| log10(stability) ~ log10(balance_f) + divergence |
240 |
5.504447 |
NA |
NA |
NA |
NA |
| log10(stability) ~ divergence |
241 |
6.364040 |
-1 |
-0.8595933 |
37.47922 |
0.000 |
Overall, balance explains more of the variance in stability than divergence, and there is virtually no difference between a model containing only balance and the full model.
Interaction divergence and richness
Richness had to be transformed to numeric and to be centered to avoid collinearity with divergence
lm_rich_div <- lm(data=complete_aggr,log10(stability)~divergence*scale(as.numeric(richness)))
# check model assumptions
# check_model(lm_rich_div)
Table 5: Type III anova table of the model with divergence and richness as predictors of stability.
anova3 <- car::Anova(lm_rich_div, type = "III")
anova_tidy3 <- broom::tidy(anova3)
# Display the tidy ANOVA table using gt with formatted p-values and adjusted size
anova_tidy3 %>%
gt() %>%
cols_label(
term = "Term",
sumsq = "Sum of Squares",
df = "DF",
statistic = "F Statistic",
p.value = "p-value"
) %>%
fmt_number(
columns = vars(p.value),
decimals = 3
) %>%
tab_style(
style = cell_text(weight = "bold"),
locations = cells_body(
columns = vars(p.value),
rows = p.value < 0.05
)
) %>%
tab_options(
table.width = px(800), # Adjust table width (e.g., 400px)
table.font.size = px(12), # Adjust font size (e.g., 12px)
data_row.padding = px(10) # Adjust row padding (e.g., 4px for more compact rows)
)
| Term |
Sum of Squares |
DF |
F Statistic |
p-value |
| (Intercept) |
11.033652088 |
1 |
442.55571734 |
0.000 |
| divergence |
0.807044347 |
1 |
32.37025122 |
0.000 |
| scale(as.numeric(richness)) |
0.001236238 |
1 |
0.04958505 |
0.824 |
| divergence:scale(as.numeric(richness)) |
0.249582101 |
1 |
10.01064605 |
0.002 |
| Residuals |
5.958668583 |
239 |
NA |
NA |
Divergence significantly interact with richness, suggesting that the relationship between divergence and stability changes with richness.
While an ideal metric of response diversity should be independent of richness.
We repeat the same model using balance instead of divergence.
lm_rich_balance <- lm(data=complete_aggr,log10(stability)~log10(balance_f)*scale(as.numeric(richness)))
# check model assumptions
# check_model(lm_rich_balance)
Table 6: Type III anova table of the model with balance and richness as predictors of stability.
anova4 <- car::Anova(lm_rich_balance, type = "III")
anova_tidy4 <- broom::tidy(anova4)
# Display the tidy ANOVA table using gt with formatted p-values and adjusted size
anova_tidy4 %>%
gt() %>%
cols_label(
term = "Term",
sumsq = "Sum of Squares",
df = "DF",
statistic = "F Statistic",
p.value = "p-value"
) %>%
fmt_number(
columns = vars(p.value),
decimals = 3
) %>%
tab_style(
style = cell_text(weight = "bold"),
locations = cells_body(
columns = vars(p.value),
rows = p.value < 0.05
)
) %>%
tab_options(
table.width = px(800), # Adjust table width (e.g., 400px)
table.font.size = px(12), # Adjust font size (e.g., 12px)
data_row.padding = px(10) # Adjust row padding (e.g., 4px for more compact rows)
)
| Term |
Sum of Squares |
DF |
F Statistic |
p-value |
| (Intercept) |
9.54928736 |
1 |
414.779284 |
0.000 |
| log10(balance_f) |
1.26534818 |
1 |
54.961191 |
0.000 |
| scale(as.numeric(richness)) |
0.02471247 |
1 |
1.073401 |
0.301 |
| log10(balance_f):scale(as.numeric(richness)) |
0.04049552 |
1 |
1.758948 |
0.186 |
| Residuals |
5.50239554 |
239 |
NA |
NA |
Balance does not significantly interact with richness, suggesting that the relationship between balance and stability is stable across richness levels.
Variable importance
Finally, we assess variable importance using the relative importance of predictors in the full model.
We use the package vip (https://cran.r-project.org/web/packages/vip/vignettes/vip.html) to calculate the relative importance of predictors in the full model.
The function vip::vip for multiple linear regression, or linear models (LMs), uses the absolute value of the -statistic as a measure of VI.
Motivation for the use of the associated 𝑡-statistic is given in Bring (1994) [https://www.tandfonline.com/doi/abs/10.1080/00031305.1994.10476059].
vip::vip(lm_div_balance)
Figure 5: Variable importance in the model including both balance and divergence as predictors of stability.
We believe that the extensive evidence here provided justifies focusing the analysis around balance, and not divergence, as a metric of response diversity.
We will thus only look at balance for the rest of the analysis.
Effect RD
We are now going to look at how response diversity (balance) affected temporal stability of total community biomass. We are going to look at the relationship between fundamental balance (so based only on species response surfaces measured in monoculture), an realised balance (measured accounting for species contribution to balance).
This is fundamentally testing our most important hypothesis.
Figure 6: Effects of fundamental and realised response diversity (measured as balance) on total community biomass temporal stability.
We can see that balance is always negatively related to temporal stability, which means that response diversity promotes stability across richness levels. Interestingly, we see that there is little difference between fundamental and realised balance. Yet, as the richness increases, the relationship between realised balance and stability becomes steeper compared to fundamental balance.
But is the difference in the slope of fundamental and realised balance significant? We can test this using a linear model with interaction between balance and type (factor: fundamental or realised balance).
Balance: realised vs fundamental
# compare if the slope of fundamental and realised balance is significantly different for each richness level
# Fit the linear model with interaction
model_F_R <- lm(log10(1/CV) ~ log10(balance) * type, data = main_r_dd)
Table 7: Type III anova table of the model with interaction between balance and type (fundamental vs realised) as predictors of stability.
| ANOVA Table for Linear Model |
| Term |
Sum of Squares |
DF |
F Statistic |
p-value |
| (Intercept) |
7.76166605 |
1 |
325.480366 |
0.000 |
| log10(balance) |
1.17164323 |
1 |
49.132089 |
0.000 |
| type |
0.09227462 |
1 |
3.869476 |
0.050 |
| log10(balance):type |
0.04940699 |
1 |
2.071850 |
0.151 |
| Residuals |
11.49415887 |
482 |
NA |
NA |
No significant difference in the slope of fundamental and realised balance was found.
Balance: realised vs fundamental by richness level
Table 7: Linear model results for the interaction between balance and type (fundamental vs realised) as predictors of stability for richness level.
| Regression Results for log10(balance):type (Realised) |
| Richness |
Term |
Estimate |
Std_Error |
T_Value |
P_Value |
| 2 |
log10(balance):typerealised |
−0.007 |
0.058 |
−0.116 |
0.908 |
| 3 |
log10(balance):typerealised |
−0.030 |
0.029 |
−1.042 |
0.299 |
| 4 |
log10(balance):typerealised |
−0.064 |
0.036 |
−1.785 |
0.076 |
Even within each richness level, the slope of fundamental and realised balance is never significantly different.
Linear models
Model: Fundamental balance
First we analyze the effect of fundamental balance, temperature, nutrients and richness on biomass temporal stability using a linear model.
balance was modelled as continuous variables, while richness, temperature and nutrients were modelled as categorical variables. balance and stability were log-transformed to meet the assumptions of linear models.
lm_full<-lm(data=complete_aggr,log10(stability)~log10(balance_f)+(richness)+nutrients+temperature)
# check model assumptions
# check_model(lm_full)
Table 8: Linear model results for the effects of balance, richness, nutrients, and temperature on community stability. Estimates are presented with 95% confidence intervals and p-values.
| Predictor |
Linear Regression Results
|
95% CI |
Linear Regression Results
|
| Estimate |
p-value |
| log10(balance_f) |
-0.05 |
-0.08, -0.02 |
<0.001 |
| richness |
|
|
|
| richness3 - richness2 |
-0.04 |
-0.09, 0.00 |
0.065 |
| richness4 - richness2 |
-0.01 |
-0.06, 0.03 |
0.8 |
| richness4 - richness3 |
0.03 |
-0.01, 0.07 |
0.3 |
| nutrients |
|
|
|
| (0.35 g/L) - (0.01 g/L) |
0.18 |
0.14, 0.22 |
<0.001 |
| (0.75 g/L) - (0.01 g/L) |
0.21 |
0.17, 0.26 |
<0.001 |
| (0.75 g/L) - (0.35 g/L) |
0.03 |
-0.01, 0.08 |
0.2 |
| temperature |
|
|
|
| (22-25 °C) - (18-21 °C) |
-0.08 |
-0.12, -0.03 |
<0.001 |
| (25-28 °C) - (18-21 °C) |
-0.10 |
-0.16, -0.04 |
<0.001 |
| (25-28 °C) - (22-25 °C) |
-0.02 |
-0.08, 0.04 |
0.7 |
A linear model was fitted to examine the effects of resource balance, richness, nutrients, and temperature on community stability (measured as log₁₀(stability)).
Among the predictors, log₁₀(balance) showed a significant negative effect on stability (Estimate = -0.05, SE = 0.016, p< 0.001). This suggests that as balance increases (more balance), stability tends to decrease.
Richness did not have a significant effect on stability within the conditions of this study.
Nutrient concentration also had a significant positive effect on stability, with estimates for 0.35 g/L (Estimate = 0.18, SE = 0.019, p < 0.001) and 0.75 g/L (Estimate = 0.21, SE = 0.019, p < 0.001) indicating increased stability with higher nutrient levels, when compared to the baseline (0.01 g/L).
Finally, temperature regimes showed a significant effect on stability. Both 22–25 °C (Estimate = -0.08, SE = 0.019, p < 0.001) and 25–28 °C (Estimate = -0.10, SE = 0.02, p < 0.001) significantly reduced stability when compared to the baseline (18–21 °C).
In summary, our findings show that temporal stability is significantly influenced by response diversity (balance), nutrient concentration, and temperature, with higher nutrient concentrations enhancing stability and higher temperatures reducing it. However, species richness was not a significant determinant of stability within the conditions of this study.
Table 9: Type II anova table of the model with balance, richness, nutrients, and temperature as predictors of stability.
| Type III ANOVA Table for Linear Model |
| Term |
Sum of Squares |
DF |
F Statistic |
p-value |
| (Intercept) |
2.11729127 |
1 |
152.932234 |
0.000 |
| log10(balance_f) |
0.15572514 |
1 |
11.248048 |
0.001 |
| richness |
0.07402585 |
2 |
2.673449 |
0.071 |
| nutrients |
1.98098126 |
2 |
71.543272 |
0.000 |
| temperature |
0.32820230 |
2 |
11.853048 |
0.000 |
| Residuals |
3.25348971 |
235 |
NA |
NA |
Interaction between temperature and nutrients
We may expect and interactive effect of the environmental variables on stability. We thus build a linear model with interaction between temperature and nutrients.
However, there is high collinearity between temperature and nutrients, which may affect the model results.
lm_full_int<-lm(data=complete_aggr,log10(stability)~log10(balance_f)+(richness)+nutrients*temperature)
# check model assumptions
check_model(lm_full_int)
So we transformed nutrients and temperature to numeric, and transformed temperature regimes in values = 1, 2, 3. Then, we centered the variables to avoid collinearity with the interaction term.
# transform nutrients and temperature to numeric. For this the units need to be removed, and temperature regimes should be transformed in values = 1, 2, 3
complete_aggr_2<- complete_aggr %>%
# Remove the units from the 'nutrients' and 'temperature' columns
mutate(
nutrients = as.numeric(gsub(" g/L", "", nutrients)), # Convert nutrients to numeric
temperature = gsub(" °C", "", temperature) # Remove the unit but keep as character
) %>%
# Convert temperature ranges to numeric codes using case_when
mutate(
temperature = case_when(
temperature == "18-21" ~ 1,
temperature == "22-25" ~ 2,
temperature == "25-28" ~ 3,
TRUE ~ NA_real_ # Handle unexpected values with NA
)
)
# Fit the linear model with interaction
lm_full_int<-lm(data=complete_aggr_2,log10(stability)~log10(balance_f)+(richness)+scale(nutrients)*scale(temperature))
# check model assumptions
check_model(lm_full_int)
Table 10: Type III anova table of the model with interaction between temperature and nutrients as predictors of stability.
| Type III ANOVA Table for Linear Model |
| Term |
Sum of Squares |
DF |
F Statistic |
p-value |
| (Intercept) |
2.07662552 |
1 |
135.208125 |
0.000 |
| log10(balance_f) |
0.06339546 |
1 |
4.127649 |
0.043 |
| richness |
0.07254258 |
2 |
2.361607 |
0.096 |
| scale(nutrients) |
1.71980335 |
1 |
111.975599 |
0.000 |
| scale(temperature) |
0.32873431 |
1 |
21.403738 |
0.000 |
| scale(nutrients):scale(temperature) |
0.02595802 |
1 |
1.690115 |
0.195 |
| Residuals |
3.62466104 |
236 |
NA |
NA |
No interaction between nutrients and temperature was found to significantly affect stability. We thus retain the simplest model without interaction term.
Asynchrony
Response diversity (aka balance) has been suggested as a mechanism that promotes temporal stability of community biomass by promoting species asynchrony.
We thus calculated the asynchrony index suggested by Gross et al. (2014)[https://www.journals.uchicago.edu/doi/epdf/10.1086/673915] to calculate the effect of asynchrony on temporal stability and to see how reponse diversity relate to asynchrony.
The index ranges between -1 and 1, with -1 indicating perfect asyncrony and 1 being perfectly synchronous, and 0 indicating random variation.
Plot stability vs. Asynchrony Gross

Figure 8: Relationship between temporal stability and asynchrony (Gross) divided by nutrient level.
The Pearson’s correlation between asynchrony and stability is significant (estimate = -0.23, p < 0.001).
cor.test((-1*async_aggr$synchrony_Gross),async_aggr$stability)
##
## Pearson's product-moment correlation
##
## data: (-1 * async_aggr$synchrony_Gross) and async_aggr$stability
## t = 3.7927, df = 239, p-value = 0.0001888
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.1153693 0.3539711
## sample estimates:
## cor
## 0.2382622
Plot Asynchrony Gross vs fundamental balance
Figure 9: Relationship between asynchrony (Gross) and fundamental balance divided by nutrient level.
The Pearson’s correlation between asynchrony and balance is significant (estimate = 18, p = 0.003).
cor.test((-1*async_aggr$synchrony_Gross),(async_aggr$balance_f))
##
## Pearson's product-moment correlation
##
## data: (-1 * async_aggr$synchrony_Gross) and (async_aggr$balance_f)
## t = -2.9796, df = 239, p-value = 0.003184
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.3082462 -0.0644258
## sample estimates:
## cor
## -0.1892515
Population stability
The relationship between community stability and the stability of the individual populations that make up the community is a key question in community ecology. Importantly, community stability can result from low population stability, if populations fluctuate asynchronously, or from high population stability, if populations do not fluctuate much.
Synthesis of the literature suggests diversity can have a positive or negantive effect on population stability (Campbell et al 2010)[https://nsojournals.onlinelibrary.wiley.com/doi/full/10.1111/j.1600-0706.2010.18768.x] and (Xu et al 2021)[https://onlinelibrary.wiley.com/doi/full/10.1111/ele.13777].
Theoretical work has suggested that community stability is a product of two quantities: the (a)synchrony of population fluctuations, and an average species-level population stability that is weighted by relative abundance (Thibaut & Connolly 2013)[https://onlinelibrary.wiley.com/doi/full/10.1111/ele.12019].
Critically, a balance value close to zero can result from high response diversity, but also from high population stability (population biomass does not change largely over time).
We want to look now at whether our new metric of balance can capture these two stabilising mechanisms.
Thus, we first calculate species-level population stability weighted by relative abundance.

–>
SEM
Finally, we use a structural equation model (SEM) to explore how stability is influenced by asynchrony, population stability, balance and, nutrient levels.
In order to develop a hypothesis regarding the influence of stability, we have drawn on existing literature. This has enabled us to posit that stability is influenced by two key factors: asynchrony and population stability. In turn, these are influenced by balance and, in our particular case, by nutrient levels.
## lavaan 0.6-19 ended normally after 1 iteration
##
## Estimator ML
## Optimization method NLMINB
## Number of model parameters 12
##
## Used Total
## Number of observations 220 241
##
## Model Test User Model:
## Standard Scaled
## Test Statistic 1.511 1.433
## Degrees of freedom 3 3
## P-value (Chi-square) 0.680 0.698
## Scaling correction factor 1.055
## Satorra-Bentler correction
##
## Model Test Baseline Model:
##
## Test statistic 628.549 624.481
## Degrees of freedom 9 9
## P-value 0.000 0.000
## Scaling correction factor 1.007
##
## User Model versus Baseline Model:
##
## Comparative Fit Index (CFI) 1.000 1.000
## Tucker-Lewis Index (TLI) 1.007 1.008
##
## Robust Comparative Fit Index (CFI) 1.000
## Robust Tucker-Lewis Index (TLI) 1.008
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) 435.602 435.602
## Loglikelihood unrestricted model (H1) NA NA
##
## Akaike (AIC) -847.204 -847.204
## Bayesian (BIC) -806.481 -806.481
## Sample-size adjusted Bayesian (SABIC) -844.509 -844.509
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.000 0.000
## 90 Percent confidence interval - lower 0.000 0.000
## 90 Percent confidence interval - upper 0.087 0.082
## P-value H_0: RMSEA <= 0.050 0.825 0.847
## P-value H_0: RMSEA >= 0.080 0.067 0.055
##
## Robust RMSEA 0.000
## 90 Percent confidence interval - lower 0.000
## 90 Percent confidence interval - upper 0.087
## P-value H_0: Robust RMSEA <= 0.050 0.831
## P-value H_0: Robust RMSEA >= 0.080 0.067
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.006 0.006
##
## Parameter Estimates:
##
## Standard errors Robust.sem
## Information Expected
## Information saturated (h1) model Structured
##
## Regressions:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## stability ~
## asynchrny_Grss 0.181 0.016 11.078 0.000 0.181 0.373
## pop_stability 1.013 0.034 29.470 0.000 1.013 0.924
## asynchrony_Gross ~
## log_balance_f -0.084 0.034 -2.493 0.013 -0.084 -0.176
## nutrients -0.199 0.025 -7.994 0.000 -0.199 -0.469
## pop_stability ~
## log_balance_f -0.064 0.008 -7.748 0.000 -0.064 -0.301
## nutrients 0.124 0.007 18.314 0.000 0.124 0.661
##
## Intercepts:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .stability 0.217 0.017 12.489 0.000 0.217 1.333
## .asynchrny_Grss -0.110 0.060 -1.817 0.069 -0.110 -0.327
## .pop_stability -0.682 0.016 -41.719 0.000 -0.682 -4.596
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .stability 0.005 0.001 9.841 0.000 0.005 0.188
## .asynchrny_Grss 0.088 0.012 7.616 0.000 0.088 0.781
## .pop_stability 0.009 0.001 9.030 0.000 0.009 0.397
##
## R-Square:
## Estimate
## stability 0.812
## asynchrny_Grss 0.219
## pop_stability 0.603
Model Fit Indices
The model fit indices suggest that the model fits the data well.
Chi-Square Test (User Model): The chi-square test statistic for the user model is χ 2 =1.626 (scaled = 1.465) with 3 degrees of freedom and a p-value of 0.653 (scaled = 0.690). This indicates a good fit, as the test is non-significant, suggesting no significant difference between the observed and model-implied covariance matrices.
Comparative Fit Index (CFI) and Tucker-Lewis Index (TLI): Both CFI and TLI values are 1.000, indicating an excellent model fit. Values close to or above 0.95 are generally considered good.
Root Mean Square Error of Approximation (RMSEA): The RMSEA is 0.000, with a 90% confidence interval ranging from 0 to 0.090 (scaled = 0.080). This indicates a very good fit, as RMSEA values below 0.05 are ideal, and values below 0.08 are acceptable. The p-values for the RMSEA hypothesis tests suggest strong support for a close fit (RMSEA <= 0.05) and little evidence for a poor fit (RMSEA >= 0.08).
Standardized Root Mean Square Residual (SRMR): The SRMR value is 0.017, which is also within the acceptable range (values below 0.08 are generally considered good).
Overall, the fit indices suggest that the model is an excellent fit for the data.
Regression Paths and Interpretation
Stability Regressions
Stability ~ Asynchrony_Gross (asynchrny_Grss): The standardized estimate for the effect of asynchrony on stability is 0.340 (p < 0.001), indicating a significant positive association. Higher asynchrony in species dynamics is associated with increased community stability.
Stability ~ Population Stability (pop_stability): The standardized estimate is 0.977 (p < 0.001), showing a strong positive relationship. This suggests that community stability is highly dependent on the stability of individual populations within the community.
Asynchrony_Gross Regressions
Asynchrony_Gross ~ Log10(Balance): The standardized estimate is -0.176 (p = 0.013), indicating a significant negative effect. Higher balance leads to lower asynchrony, suggesting that as balance increases, species within the community fluctuate more synchronously.
Asynchrony_Gross ~ Nutrients: The standardized estimate is -0.469 (p < 0.001), showing a strong negative relationship. Higher nutrient levels appear to reduce asynchrony, possibly by causing similar responses across species.
Population Stability Regressions
Population Stability ~ Log10(Balance): The standardized estimate is -0.296 (p < 0.001), indicating that higher balance is associated with lower population stability.
Population Stability ~ Nutrients: The standardized estimate is 0.635 (p < 0.001), showing that higher nutrient levels are associated with increased population stability, likely because nutrients enhance conditions that support stable population dynamics.
Variances and R-Squared Values
R-Squared for Stability: The model explains 90.4% of the variance in community stability, indicating strong predictive power.
R-Squared for Asynchrony_Gross: The model explains 21.9% of the variance in asynchrony, which is moderate.
R-Squared for Population Stability: The model explains 56.2% of the variance in population stability, showing that nutrients and balance are important but not the only factors influencing it.
Summary Interpretation
Model Fit: The model has an excellent fit, as indicated by the fit indices.
Stability: Community stability is strongly influenced by both population stability and asynchrony among species, with population stability being the stronger predictor.
Asynchrony and Balance: Asynchrony decreases with increasing balance and nutrients, suggesting that these factors promote more synchronized fluctuations among species.
Population Stability and Nutrients: Higher nutrient levels are associated with increased population stability, suggesting that nutrient availability supports stable population dynamics. Conversely, higher balance is associated with decreased population stability.